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by Philip R. Nachtsheim and Paul Swigert 
Lewis Research Center 

SUMMARY 

A method for the numerical solution of differential equations of the boundary-layer 
type is presented. The asymptotic boundary conditions are satisfied at the edge of the 
boundary layer by adjusting the initial conditions so that the mean square error between 
the computed variables and the asymptotic values is minimized. Convergence to a solu- 
tion is rapid and appears to be insensitive to the initial guesses of the initial conditions. 
Use of a least-squares convergence criterion leads to the unique solution even in the 
case of retarded main stream flows. A description of the method is given, and two 
examples taken from boundary -layer theory are worked. For the Falkner-Skan differen- 
tial equation, the relation between the wedge angle and the shear stress at the wall is 
obtained. 


INTRODUCTION 

The numerical integration of the ordinary differential equations of boundary-layer 
theory involves the satisfaction of asymptotic boundary conditions; that is, some bound- 
ary conditions are specified at the initial point or wall, and others are specified as limits 
that must be approached at large values of the independent variable corresponding to the 
edge of the boundary layer. In order to integrate the differential equations numerically 
from the wall to the edge of the boundary layer, it is necessary to specify as many addi- 
tional conditions at the wall as there are conditions to be satisfied at the edge of the 
boundary layer. These additional initial conditions have to be varied in order to satisfy 
the conditions at the edge of the boundary layer. All methods that have been used to find 
the proper initial conditions rely on the fact that, for large values of the independent 
variable, the integrals of the differential equations depend on the initial conditions. One 
method that has been used to find the proper initial conditions, called the initial-value 
method (ref. 1), is that of obtaining integrals of the differential equations with guessed 
initial conditions; that is, an attempt is made to integrate the differential equations to 



the edge of the boundary layer. If this can be accomplished, corrections are made for 
the initial guesses by the Newton-Raphson method, and the process is repeated until con- 
vergence is achieved. 

In order to carry out this iteration, additional differential equations have to be inte- 
grated. The additional differential equations are obtained by differentiating the terms of 
the original differential equation with respect to the initial conditions. The integrals of 
these additional differential equations, henceforth referred to as perturbation equations, 
give the rate of change of the integrals of the original differential equation with respect 
to the initial conditions. 

Another method that has been used, called quasilinearization (ref. 2), is similar in 
principle to the initial-value method, except that the original differential equations are 
linearized. The resulting linear differential equations for the current approximation are 
inhomogeneous, since they also contain the previous approximation as members. 

The integrals of the perturbation equations of the initial-value method correspond to 
a complementary solution of the linear differential equations of the quasilinearization 
method. In application, the two methods are different in that, in the quasilinearization 
method, the complementary solutions are used to obtain solutions to the inhomogeneous 
linear differential equations, whereas in the initial-value method, the integrals of the 
perturbation equations are used to adjust the starting values at the initial point. A dis- 
advantage associated with the quasilinearization method is that functions representing the 
previous approximation have to be stored in order to construct the current approximation. 
Furthermore, the solution of the inhomogeneous equation obtained by combining com- 
plementary solutions and a particular integral is usually not well determined except near 
the origin or initial point, as pointed out by Hartree (ref. 3). The solutions of the per- 
turbation equations, however, can be used to determine the proper initial conditions 
closely. This is the procedure that is followed in the initial-value method. 

Both methods have been used to obtain solutions of the boundary-layer equations. 
Neither of these methods., however, has solved the problem of when to stop the integra- 
tion. Since one boundary condition is specified as a limit that must be approached at 
large values of the independent variable, the integration should be stopped at a value of 
the independent variable that is sufficiently large so that the various functions are 
approaching their asymptotic values. This value of the independent variable could 
properly be called the edge of the boundary layer. In the various methods used to solve 
this problem, this value of the independent variable must be guessed beforehand. If the 
guessed value is too small, the integrals of the differential equation may not be able to 
satisfy the imposed conditions; or, there may be more than one value of an initial condi- 
tion that leads to integrals of the differential equations that satisfy the imposed boundary 
conditions. 

If the guessed value is too large, it is possible that the integrals of the equations 
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will diverge, or that convergence to a solution is extremely slow. 

Another problem that sometimes besets the numerical integration of the boundary- 
layer equations is the apparent insensitivity of the integrals of the boundary-layer equa- 
tions to the initial conditions. This difficulty appears in the integration of the Falkner- 
Skan differential equation (ref. 4). All the integrals of the Falkner-Skan equation for the 
case of nonretarded flows in the main stream tend to diverge except the ones with the 
proper initial conditions; however, for the case of retarded flows, a solution with any 
value of the initial condition near the correct one will eventually meet the proper boundary 
condition at infinity. The unique solution was found by Hartree (ref. 4) for retarded flows 
by imposing the additional condition that the desired solution is the one that approaches 
the boundary condition most rapidly from below. However, this behavior of the solutions 
can only be learned from several trial runs. 

In an effort to adapt the initial-value method to the solution of differential equations 
with asymptotic boundary conditions, a method was developed to eliminate the problem 
of when to stop the integration. This method of solution is capable of satisfying the 
boundary conditions at the edge of the boundary layer correctly; that is, the boundary 
values are approached asymptotically. This is accomplished by choosing the additional 
initial conditions so that the mean square error between the computed variables and the 
asymptotic values is minimized. A method of solution is proposed that is based on a 
least-squares convergence criterion in which the edge of the boundary layer is approached 
in steps. The convergence rate to a solution is rapid, and convergence appears to be 
insensitive to the initial guesses of the initial conditions. Use of a least-squares con- 
vergence criterion leads to the unique solution even in the case of retarded main stream 
flows. A description of the method is given, and two examples taken from boundary- 
layer theory are worked. By using this method of solution it is feasible to obtain directly, 
for those differential equations which contain a parameter, the relation between the 
properties of the solutions and the values of the parameter. For the Falkner-Skan differ- 
ential equation, the relation between the wedge angle and the shear stress at the wall is 
obtained. 


DESCRIPTION OF METHOD 

A description of the method can best be given by referring to a definite example of a 
boundary -value problem with an asymptotic boundary condition. In Falkner and Skan's 
treatment of the laminar boundary layer of an incompressible fluid, the following equa- 
tion arises (see ref. 4): 


f ,M = -ff" + |6(f' 2 - 1) 


( 1 ) 


3 



with the boundary conditions 


77 = 0: f = f* = 0 

77 — f' — 1 

(Symbols are defined in appendix A. ) Primes denote differentiation with respect to 77 , a 
measure of the distance from the wall. The dependent variable f is related to the usual 
stream function. The potential flow is given by a power law. The flow velocity outside 
the boundary layer in the main stream is proportional to the distance along the wall 
raised to the power 0/(2 - 0 ). For retarded main stream flow 0 < 0, and for non- 
retarded flow 0 > 0. 

In practice, the asymptotic boundary condition is replaced by the condition that 
f’ = 1 to a sufficient degree of accuracy at 7 ? = 77 edg e , where ? 7 edge is the value of the 
independent variable at the edge of the boundary layer. The boundary -value problem is 
equivalent to the problem of finding a value of f M (0) for which the boundary condition at 
the edge of the boundary layer is satisfied; that is, it is desired to find a solution of the 
nonlinear equation 


= 1 

where f^ dge = f'(?? edge )- The function P dge [f"(°)] in this problem is not, in general, 
explicit; it will be expressed as a function of f M (0) through an integration of equation (1). 
With the notation x = f”(0), observe that a small change Ax in x changes f' by the 
amount 


3x 

so that the necessary correction to a first approximation comes from the solution of the 
linear equation 


1 = f T + — Ax 
3x 


at * = ^edge' 

The solution of the equation for Ax can be performed provided that the partial 
derivative of f' with respect to x can be evaluated at ?7 ed g e - The partial derivative 
can be evaluated by forming the perturbation equation for the function f'. This equation 
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f'(0) Distance from wall, tj 


(a) r'(0) at 7) - 5. (c) f' for r'<0) - 0. 85. 

Figure 1. - Variation of f and f". Falkner-Skan parameter, L 

is obtained by differentiating the terms in equation (1) appropriately. With the notation 



0f ft _ Sf' 

1 > • • • 

dx x 3x 


the following perturbation equation is obtained: 

f x’ = " (ff x + f ' ,f x) + 2/3f ’ f x ( 2 ) 

with the initial conditions 


0: f x = f x=°’ f x = 1 


Given an initial estimate of x = f’’(0), subsequent values of x can be computed by 
integrating equations (l)and (2) and applying the Newton-Raphson method to obtain correc- 
tions. 


The problem of where to stop the integration remains; that is, r 7 ec jg e i s unknown. 
Figure 1(a) illustrates the difficulties that arise when the boundary condition f' = 1 is 
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applied at a finite value of rj. The data for figure 1(a) were obtained by integrating 
equation (1) for specified values of f"(0) with /3 = 1. The integration was stopped at 
rj = 5. Tentative values of ?? ec jg e hereinafter be referred to as V s t 0 p It should be 
noted that the problem is not that V e ^g e i- s unknown, but an appropriate V s ^ 0 ^ must be 
determined in some fashion for use in the as yet unconverged iterations. 

Figure 1(a) shows f' and f" as a function of x= f"(0) at V s t Q p = 5. It can be 
seen that there are two values of x that satisfy the condition f' = 1 when 7? st0 p is 
taken to be 5, namely, approximately x = 0. 85 and x = 1. 23. Figures 1(b) and (c) show 
the two solutions corresponding to the two values of x. It can be concluded from this 
illustration that the satisfaction of the single boundary condition f' = 1 at 77^ = 5 

does not lead to a unique solution. The proper value of x can be determined by ob- 
serving where f" = 0. From figure 1(a) it can be seen that the correct value of x is 
approximately 1. 23. Satisfaction of both boundary conditions f' - 1 and f" = 0 at 
?7stop = 5 ensures the satisfaction of the original boundary condition of f’, namely, that 
f* approaches unity asymptotically as r) approaches infinity. Hence, in order to satisfy 
the asymptotic condition at a finite value of t], the conditions f' = 1 and f" = 0 should 
both be imposed (examination of the differential equation, (eq. (1)), will reveal that all 
higher derivatives will be zero also); that is, Ax should be chosen so that both equations 


1 = f» + f’ Ax 

X 


(3) 


and 


0 = f" + f” Ax 


(4) 


are satisfied at 77 = h s t 0 p- This is, in general, impossible since there are two condi- 
tions and only one adjustable parameter Ax. However, a satisfactory solution that is 
consistent with the idea that the boundary condition cannot be satisfied exactly at a finite 
value of 7 ] is to seek the least-squares solution of equations (3) and (4) (ref. 5). Hence 
consider the discrepancies 6^ and fig, where 


f’ Ax + f' - 1 

X 


5 9 = f' ' Ax + f" 

2 2 

and attempt to find a value for Ax that will make the sum 6^ + 62 as small as possible. 
To minimize the sum, it is necessary to equate its first derivative with respect to Ax 
to zero. This calculation yields a value of Ax corresponding to the minimum of the sum 
Sj + Sg, that is, 
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Figure 2. - Error as function of f 1 10) for different values of T) st0 p. Falkner-Skan param- 
eter, 1. 


Ax = 


f^U - f') - f”f" 




( 5 ) 


Along with the discrepancies 6 ^ and 62 consider the sum of the squares of the 
deviations of the computed quantities from their asymptotic values 

E = (1 - f ’) 2 + f ” 2 ( 6 ) 


Henceforth, E shall be referred to as the error. The magnitude of E evaluated at 
^stop &i ves an indication of how unsatisfactory the asymptotic boundary conditions are. 
The value of x that gives Ax = 0 corresponds to the minimum with respect to x of E, 
as can be seen from equation (5). In figure 2, the quantity E is shown plotted as a 
function of x for various values of ?7 st 0 p- The data for figure 2 were also obtained 
from the integrations that were used to prepare figure 1. It can be seen that the least- 
squares solution gives a unique solution for x. It is this property, that the minimum of 
E with respect to x, E min , corresponds to no change in the required initial conditions, 
that permits the asymptotic boundary condition to be approximately satisfied at a finite 
value of the independent variable. 

For this example, the Falkner-Skan solution with = 1, it appears that the quantity 
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E has a definite minimum at each value of q, even at q = 0, as seen from figure 2. The 
edge of the boundary layer could be defined as that value of q for which the value of 
E m in i s * ess than some preassigned small value as the range of integration is increased. 
This gives a value of V e ^g e - 

The concept of convergence in the least-squares sense allows the initial-value 
method to be modified in order to solve a wide class of boundary -layer problems. To 
understand why the initial-value method requires modification, consider how the method 
can fail when it is applied to boundary-layer-type problems. 

Recall that in the initial-value method corrections are obtained to the guessed 
initial conditions after integrating the differential equations to the edge of the boundary 
layer. This method can fail if the initial guess is so poor that the integration "blows up" 
(i. e. , diverges beyond a prescribed limit) before reaching the edge of the boundary 
layer. 

This problem can be avoided and corrections to the initial conditions can be obtained 
by attempting to minimize the mean square error between the computed solution and the 
asymptotic values before reaching the edge of the boundary layer. The modification of 
the initial-value method consists then of finding E m ^ n at any arbitrary value of 77 ^ . 

By this modified procedure, the choice of ? 7 st 0 p, where the minimizing process is first 
carried out, may be dictated by other considerations than where the edge of the boundary 
layer is, such as keeping the variables within prescribed limits so that the solution does 
not diverge. After E m ^ n has been found, no further corrections to the initial conditions 
need to be made at the chosen value of ?? st 0 p, and the range of integration can be in- 
creased if necessary. 

If the value of E m ^ n found is not small enough, it will be necessary to increase the 
range of integration. In this way, the edge of the boundary layer is approached in steps 
by successively increasing the range of integration until the correct initial conditions and 
the value of ? 7 ed g e are found. In other words, as f? st0 p is increased, E min tends to 
zero, and the required initial condition x = f"(0) also tends to a definite value. This 
gives the solution of the boundary-value problem. 


NUMERICAL SOLUTION 

The method described in the previous section was programed in FORTRAN IV, 
double precision on the IBM 709411-7040 direct-couple system. The boundary-layer and 
perturbation differential equations were rewritten as systems of first-order differential 
equations and integrated with a predictor-corrector (Adams -Moulton) subroutine using 
one correction per step and a fixed increment. 

The boundary-layer equations along with the perturbation equations were integrated 
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Figure 3. - Convergence history of solution of equation (1) 
Falkner-Skan parameter, 1. 


to the specified value of the independent vari- 
able. The corrections were then determined, 
and the process was repeated until the relative 
change in the correction term, equation (5), 
was less than a small preassigned value. 

Extreme accuracy was not required at the 
smaller values of ?7 st0 p since the boundary 

conditions cannot be well satisfied there. A 
-8 

value of 1X10 was used at the larger values of 
77 g top check the relative change in the correc- 
tion term. This small value seemed reasonable 
since convergence is so rapid. 


RESULTS OF NUMERICAL SOLUTION FOR EXAMPLE 1 


Figure 3 shows the solution f' of equation (1) with ft = 1 as a function of the 
independent variable t], for the five trials needed for convergence. With an initial guess 
of 1. 0 for f”(0), the equations were integrated to a value of V s i Q p = 2 twice and to a 
value of 7? s j. 0 p - 5 three times. Observe that, for the first integration, the values of f' 
tend to deviate radically from the correct solution at a relatively small value of rj. If 
these calculations had been carried out to a larger value of 77, the values of f’ would 
have attained a magnitude so large that any use of these large numbers in a Newton- 
Raphson scheme would be meaningless. 


TABLE I. - EFFECT OF ?j gtop ON CONVERGENCE RATE 


f”(0) 

^stop 

2 

^stop 

5 

(initial guess) 






f”(0) 

Percent 

f”(0) 

Per cent 


(after two trials) 

difference 

(after two trials) 

difference 

0.25 

1.2799034 

3.83 

0.26706180 

-78.33 

. 50 

1.2449846 

1.00 

. 69205364 

-43.85 

.75 

1.2292604 

-.26 

. 88486032 

-28.21 

1.00 

1.2266764 

-.47 

1.0413712 

-15.51 

1.25 

1.2266282 

-.48 

1.2325888 

.00 

1. 50 

1.2266765 

-.47 

1. 2444813 

.96 

1.75 

1.2277231 

-.39 

1. 2890732 

4.58 

2.00 

1.2312423 

-. 10 

1.3564389 

10. 04 

2.25 

1.2383960 

.47 

1. 8834930 

52.80 

2. 50 

1.2498889 

1. 40 

2. 4946776 

102.39 

2. 75 

1 2660459 

2 71 



3.00 

1.2869264 

4. 40 
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TABLE H. - CONVERGENCE HISTORY OF INITIAL VALUES 


(a) Example 1 


Trial 

Correction 

^stop 

f"(0) 

1 


2.0 

1.0 

2 

1 

2.0 

1. 2463981 

3 

2 

5.0 

1.2266764 

4 

3 

5.0 

1. 2326729 

5 

4 

5.0 

1. 2325878 


5 


1.2325878 


Literature values 


Source 


Ref. 2 
Ref. 6 


^stop 

5.0 

4.0 


f"(0) 

1. 2397 
1. 2325876 


(b) Example 2 


The effect of V s ^- 0 p on the convergence 
process is illustrated in figure 2 and table I. 
Figure 2 shows the error, equation (6), as 
a function of f"(0) for different values of 
^stop" Observe that by initially using a 
small value for "stop (’'stop = 2 in 2) 
the range of meaningful initial guesses is 
increased, that is, initial guesses that 
yield a relatively small error. For ex- 
tremely large values of V s i Q p the 
parabola-like curves in figure 2 degenerate 
to a vertical line, and the initial guess is 
limited to the correct answer. 

Table I shows the effect of ?7 st0 p on 
the convergence rate. Also, the values of 
f”(0) after two trials and the percentage 
difference from the correct value of f"(0) 
for a wide range of initial guesses and two 
values of V s i 0 p are given. The conver- 
gence to the correct answer is always 
faster for the smaller V s ^ 0 p except when 
the initial guess is close to the correct 
answer. In fact, the two largest initial 
guesses ’’overflowed” the computer before 
the larger 77^ was reached. This indi- 
cates that a very crude initial guess may 
be used if the integration is carried out to 
a small 77 stop . 

T' give an idea of the rate of conver- 
gence to a solution for this example using an initial guess of f”(0) = 1, table n(a) is pre- 
sented. Using quasilinearization and five corrections, Radbill (ref. 2) obtained a value 
of f’’(0) that agrees with that of Yohner and Hansen (ref. 6) to within seven units in the 
third decimal place. The method developed herein reached an answer after two correc- 
tions that agrees with that of Yohner and Hansen to within six units in the third decimal 
place. After two more corrections the answer agreed to within two units in the seventh 

decimal place (see table 11(a)). The last integration was performed as a stopping criterion 

—4 

for the computer. Computer time for the complete solution, using a step size of 2 was 
approximately 0. 03 minute. 

This example was also programmed by single precision with good results. 


Trial 

Correction 

^stop 

f"(0) 

h'(0) 

1 


2.0 

1 . 0 

- 1 . 0 

2 

1 

2.0 

.61735605 

-. 59008297 

3 

2 

8.0 

.63211384 

-. 57655550 

4 

3 

8.0 

. 66795925 

-.52593128 

5 

4 

8.0 

. 67387305 

-.50916398 

6 

5 

8.0 

. 67412049 

-.50790411 

7 

6 

8.0 

. 67412438 

-.50789273 


7 


. 67412438 

-. 50789273 


Literature values 


1 

Source | 

i 

" -op 

f’’(0) 

h'(0) 

Ref. 8 


8.0 

0.6741 

-0. 5080 


The 
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difficulty in using single precision presented itselt in the form of reducing the range of 
initial guesses that yielded a relatively small error. Evidently, this was caused by the 
increased round-off error associated with the use of single precision. 

Equation (1) was also integrated for the case of a retarded main stream flow with 
/3 = -0. 1. However, no additional conditions such as those used by Hartree (ref. 4) had to 
be imposed in order to obtain the unique solution. The value of f”(0), which was obtained 
by minimizing the mean square error, agreed with the value obtained by Smith (ref. 7) 
for this case to six decimal places. A description of the flow chart and the FORTRAN 
listing for this program are presented in appendix B. 

EXTENSION TO SYSTEMS OF EQUATIONS 

The least-squares modification of the initial-value method described in the preceding 
section for a differential equation with a single dependent variable can easily be general- 
ized to systems of equations. As an example of a problem with two dependent variables, 
consider the boundary-value problem for the free-convection flow about a vertical plate. 
This example is taken from reference 8. It consists of solving the differential equations 

f,T> = + 2f' 2 - h (7) 

h” = -SPrfh’ (8) 

with the boundary conditions 

n = 0: f = f» = 0, h = 1 
77 — °°: f’ -«• 0, h — 0 

Again the primes denote differentiation with respect to 77 (a measure of the distance 
from the wall), and the dependent variable f is a dimensionless form of the usual stream 
function. The dependent variable h is proportional to the temperature excess of the 
fluid over the ambient temperature, and Pr denotes the Prandtl number. 

Since there are two asymptotic boundary conditions to satisfy, two additional initial 
conditions at the wall will have to be adjusted; that is, values of x s f"(0) and y = h’(0) 
are sought that will satisfy simultaneously the nonlinear equations. 

h edge[ f,, (°)’ h, (°)] = 0 

and 

f edge[ f,,( °)’ = 0 

where f^ dge = fledge* and h edge s ^edge** The necessar y first corrections to a 
first approximation x and y come from a solution of the linear equations 
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and 


0 = f' + Ax + Ay 


0 = h + hxAx + hyAy 

at rj = J? edge , where f^ = df'/dx, = 3f'/9y, etc. However, in order to satisfy the 
asymptotic boundary conditions, the preceding equations must be supplemented by 


and 


0 = f" + Ax + f" 


Ay 


0 = h' + h^ Ax + h^ Ay 

and Ax and Ay must be found such that the sum of squares of the errors in the pre- 
ceding four equations be a minimum. The least-squares solution for those quantities is 
given by the solution of the following matrix equation: 


f x 2 +h l +f x 2 + h x 2 ’ 

Vy + W + f i' f y' 

_ f x f y + h x h y + f, x f ’y + h x h y> 

f y 2 + *V + f y 2 + ^ 




Ax 


Ay 


f’f' + hh + f"f" + h'h' 

XXX X 

f*r + hhy + f”fy + h’hy 


O) 


The values of x and y that give Ax = Ay = 0 correspond to the minimum with 
respect to x and y of the quantity 

9 9 9 9 

E = r + h + f”^ + tr (10) 

The partial derivatives with respect to x and y that appear in equation (9) are ob- 
tained by integrating the appropriate perturbation differential equations. The perturba- 
tion differential equations for the x derivatives are 

f x' = _3(f x fM + ff x> + 4ff x “ h x Ul) 

and 

h” = -3Pr(f x h’ + fh^) (12) 

with the initial conditions 


r? - 0: f = f* = h = h' = 0: f" = 1 

1 X X X X ’ X 
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' 0 1 2 3 4 5 6 7 8 

Distance from wall, rj 

Figure 4. - Convergence history of solutions of equations (7) and (8). Prandtl number, 0. 733. 

The perturbation differential equations for the y derivatives are 

f”' = -3(f y f" + ff”) + 4f’f^ - hy (13) 

and 

h^' = -3Pr(f y h’ + fhp (14) 

with the initial conditions 

^0: fy = f , =f ,, = hy = 0; h . =1 

Note that the equations for the y derivatives are the same as the equations for the 
x derivatives. The integrals of these equations will differ since the initial conditions 
are different. 

Note also that there are as many additional systems of perturbation equations to 
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0 


-.2 
-.4 

h'(0) 6 

-.8 
- 1.0 
- 1.2 

.2 .4 .6 .8 

f'(0) 


E > 1. 0 


+ f"2 + h'2J 


1.0 1.2 .4 .6 .8 1.0 1.2 

f"(0) 


^ ^stop = 2. (b) ^ stop " 

Figures. - Level curves of error against f"(0) and h'(0). 


integrate as there are asymptotic boundary conditions to meet. 




NUMERICAL SOLUTION OF CONVECTIVE FLOW PROBLEM FOR EXAMPLE 2 


The same general procedure used in the previous example was employed to solve 
this example. One notable difference is that five trials were needed at the larger value 
of ?? s t 0 p, rather than the three trials of the first example. This is probably because of 
the need of adjusting two initial conditions rather than one and because an error in one 
initial condition leads to an error in computing the correction term of the other. 

The solutions of equations (7) and (8) are shown in figure 4 for Pr = 0. 733. Note 
(as in the previous example) that the functions f’ and h of figure 4 diverge radically 
from the true solutions, for the initial guesses of trial 1. 

Table 11(b) (p. 10) shows the rapid rate of convergence for this more difficult ex- 
ample. With initial guesses of 1. 0 and -1.0 for f”(0) and h' (0), respectively, con- 
vergence is realized after two trials with ?? sto p = 2 and five trials with T7 stop = 8. The 
results are in close agreement with the published values of Ostrach (ref. 8). 

The effect of V s i 0 -p c» n the convergence process for this example is illustrated in 
figures 5(a) and (b). Level curves for the error, equation (10), are plotted in the 
f"(0),h'(0) plane for = 2 in figure 5(a), and for ^ - 5 in figure 5(b). As in 

the previous example, by using a small value for V s t 0 p> the range of initial guesses 
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that yields a relatively small error is again increased; that is, the area where E is 
less than 1 is much greater in figure 5(a) than in figure 5(b) . For extremely large values 
of 77 s t 0 p) that is, V s i Q ^ greater than 8, the area enclosed by the level curve where E 
is less than 1 shrinks to a point, and the initial guess is limited to the correct answer. 

Computer time for the second example was about 0. 10 minute with a step size of 

— 4 

2 . In both examples, the convergence was made insensitive to the initial guesses by 

choosing a small value of V s i 0 -g for the first two trials. The answer thus obtained for 
the initial conditions was good enough so that V s t Q -g could be increased considerably. 

The error term may be made as small as desirable by continuous extension of the range 
of integration. This error term may also be used as a stopping criterion for computer 
calculations and finding V e ^g e automatically. A discussion of the flow chart and the 
FORTRAN listing for this program are given in appendix B. 


SOLUTIONS OF AN EQUATION CONTAINING A PARAMETER 

Very often it is desired to examine the general properties of the solutions of a 
differential equation containing a parameter when the parameter is varied. For example, 
in the case of the Falkner and Skan differential equation (eq. (1)), which contains the 
parameter P, it is desirable to know the variation of f”(0) with /3. The physical 
significance of f"(0) is that it is proportional to the shear stress at the wall. 

If /3 > 0, the flow may be regarded as that near a wedge, and the included angle of 
the wedge is Pit. But if P < 0, it can only be said that it corresponds to a flow in an 
adverse pressure gradient (retarded main stream flow). It is desired to know the be- 
havior of the shear stress as P is varied. It is especially important to know what value 
of P gives f"(0) = 0, that is, the laminar separation profile. This behavior of the solu- 
tions can be learned by solving equation ( 1) for a preassigned set of values of P and 
determining a set of values of f"(0). 

Since the modification of the initial-value method described previously allows accu- 
rate solutions of equation (1) to be obtained quickly and easily, it is feasible to obtain the 
variation of f"(0) with P by using a more economical procedure than by obtaining solu- 
tions for a preassigned set of values of P. This procedure is essentially the method of 
steepest descent, as explained in reference 9. It consists of finding the curve of solu- 
tions f"(0) against P in the P, f"(0) plane and following it. The points that lie on the 
curve of solutions are the values that allow equation (1) to satisfy the imposed boundary 
conditions. 

The curve of solutions of equation (1) in the P, f”(0) plane can be considered to be 
given implicitly by the equation 
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(15) 


' 1 

where f edge ~ f ’^edge^' ^ P oint /3, f"(0) in the plane may be regarded as lying on a 
curve fg^gg = constant. If /3, f M (0) is an approximate solution of equation (15), the curve 
f^dge = constant will be nearly parallel to the curve of solutions f^ge = 1, and the 
necessary corrections to a first approximation can be found by solving the equation 

1 = f’ + fg A/3 + f^ Ax (16) 


together with the equation 


0 = -f^ A/3 + fjg Ax 


(17) 


at V = 7 e dge’ where f £ = 3F/3/3 and f^ = df'/dx, x = f”(0). Equation (17) guarantees 
that the increment (A/3, Ax) will be perpendicular to the curves fg dge = constant. Succes- 
sive iterations lead to a point on the curve f^ = 1. The direction of the tangent can be 
obtained by solving 


° = f £ +f x^ ( 18 ) 

for dx/d/3 and may then be used to predict further points close to the curve. 

The solution of the equations for (A/3, Ax) (eqs. (16) and (17)) can be performed pro- 
vided that the partial derivatives of f' with respect to x and /3 can be evaluated at 
77 e dg e - The partial derivatives with respect to x can be obtained from the solution of 
equation (2), as explained previously. The partial derivatives for the /3 variation can be 
evaluated by forming the perturbation equation for the function f'. This equation is ob- 
tained by differentiating the terms in equation (1) appropriately with respect to /3. The 
following perturbation equation for the /3 variation is obtained 

f /3" = "( ff /3 + f ’V + 2 # ,f /3 + (f' 2 “ 1) (19) 


with the initial conditions 


7-0: f fi = = fjj = 0 

Note that equation (19) is an inhomogeneous differential equation so that the homogeneous 
initial conditions do not necessarily yield a trivial solution. 

Furthermore, note that, if the requirement of steepest descent is dropped, for 
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example, that is, the curve of solutions is approached on a line parallel to the £-axis, 
equation (19) alone could be used in conjunction with equation (1) without the use of equa- 
tion (2) to obtain the curve of solutions for preassigned values of x = f”(0). In this way, 
the value of /3 that gives the separation profile f”(0) = 0 can be found. 

The requirement of the asymptotic approach of the solutions to the proper values is 
incorporated by recognizing that the curve of solutions of equation (1) can also be con- 
sidered to lie on the curve given implicitly by the equation 

f edge[ /3 ’ f,,{0) ] =0 ( 2 °) 

and requiring that this equation, where Pg<jg e = f" (T? e dge^’ satis f iec ^ as we ^ as equa- 
tion (15) at r/ = ri e £ . In this case, if /3, f"(0) is an approximate solution of equa- 
tion (20), the necessary corrections to a first approximation can be found by solving the 
equation 


0 = f” + fjj A/3 + f" Ax 


( 21 ) 


together with the equation 


0 = -f” A/3 + f£ Ax (22) 

at 7 ] = ?7 e( ig e ' The direction of the tangent can be obtained by solving 

0 = f !d + f ' ’ — (23) 

^ x d/3 

for dx/d/3. The required values of A/3, Ax for the asymptotic approach are obtained 
from the least-squares solution of equations (16), (17), (21), and (22). The least-squares 
solution for these quantities is given by 


A/3 = 


fj^P - 1) + fjgP’ 

f /3 2 + f x 2 + f /3 2 + f x 2 


(24) 


and 


f^(P - 1) + f”f” 

Ax = 

l h 2 + * l t + f x 2 


( 25 ) 
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and the least-squares solution of 
equations (18) and (23) gives the direc- 
tion of the tangent 


dx 

d/3 


vi : gg 

f i 2 +i x 2 


(26) 


The process just described was 
programed for solution on the com- 
puter. Very few changes were needed 
in the program used in example 1. 
Appendix B gives a discussion and the 
FORTRAN listing of this program. 

The entire curve of solutions 

(fig. 6), defined by 23 points over the range /3 = -0. 1 to 2. 1, was obtained in approxi- 
mately 0. 8 minute of computer time. Two values of ?7 sto p were used for this solution; 
they are 4. 0 and 8. 0. The smaller r? sto p was used only twice to start the solution since 
the predicted points on the curve of solutions, obtained from equation (26), were accu- 
rate enough to integrate directly to the larger ?7 sto p* A solution was obtained in about 
four iterations, for each point, in this manner. 

Equations (2) and (19) were integrated along with equation (1) to give the necessary 
terms for the corrections obtained from equations (24) and (25). When convergence was 
realized, equation (26) was applied to predict a new solution. 

The value of /3 where f"(0) = 0 was obtained by the method of example 1 by treat- 
ing /3 as the dependent variable; that is, f"(0) was held constant, while /3 was allowed 
to vary. This is accomplished by integrating equation (19) along with equation (1) and 
applying the correction factor given by 


iUl - f') - fgf* 

A/3 = A P 


ft 2 « r 2 

P + P 


Note the similarity to equation (5). 

-28 

This method with V 8 ^ 0 = 16 yielded an error term (eq. (6)) of less than 1X10 
The numerical value of /3 when f"(0) = 0 is -0. 198837682 and agrees with that of 
reference 7 to six decimal places. When the value of /3 presented in reference 7 was 
used in the program of example 1, the solution oscillated, and convergence could not 
be achieved. When the value presented herein was used, no such oscillation was 
experienced. 
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CONCLUSIONS 


The two main problems of integrating the boundary -layer equations, namely, approxi- 
mating the missing initial conditions and determining when to stop the integration in the 
as yet unconverged iterations, have been reduced to an automatic initial-value technique 
that is easily programed on high-speed computers. The technique was successfully 
applied to two problems taken from boundary- layer theory and offers promise of producing 
equivalent results when applied to other boundary-layer problems. The method appears 
to be insensitive to initial guesses and converges quickly to the solution. By this method 
the implicit study of differential equations containing a parameter has also been reduced 
to an automatic technique. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, June 11, 1965. 
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APPENDIX A 


SYMBOLS 


E sum of squares of errors in 
boundary conditions 

f related to stream function 

h related to temperature excess of 

fluid over ambient temperature 

Pr Prandtl number 

x x = f"( 0 ) 

Ax increment in x 

y y = h’(o) 

Ay increment in y 

/3 Falkner-Skan parameter 

A/3 increment in /3 

measure of distance from wall 


Subscripts: 

edge edge of boundary layer 

stop tentative values of edge of boundary 
layer 

x denotes partial differentiation with 

respect to x 

y denotes partial differentiation with 

respect to y 

/3 denotes partial differentiation with 

respect to j3 

Superscript: 

' denotes differentiation with respect 

to 77 
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APPENDIX B 


DESCRIPTION OF FLOW CHART OF MAIN PROGRAM 


A general flow chart (fig. 7) is presented for the three main programs listed. The 
numbers in the following description refer to the flow chart box numbers and formula 
numbers of the FORTRAN listings. 

1. Input, such as first guess of initial conditions, frequency of output, increment of 
independent variable, various end points of the independent variable, and convergence 
tests, is read into storage. 

2. The title of the problem and the input may be printed at this step. 

3. Initial values of equations are set in this box. This includes first guesses also. 

4. Solution is advanced one increment in the independent variable direction. This 
is done by a subroutine that uses either Runge-Kutta or Adams- Moulton integration 
schemes. This subroutine is described in appendix C. 
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Figure 7. - General flow chart for MAIN program. 


5. Output of the dependent vari- 
ables may be written at this point. A 
test of the independent variable may be 
used so that output is not written at 
every increment. 

6. If the independent variable has 
not reached the end point, the solution 
is advanced another increment. 

7. Changes in starting values are 
computed by using equations presented 
in the text. Starting values are ad- 
justed. The error term is computed. 

8. Output of starting values is 
performed here, if desired. 

9. Changes in starting values are 
tested to determine if convergence is 
realized. If not, iteration is repeated. 

10. Test of the error term is made 


and. 


is increased if the error 


^stop 

does not pass the test. 

11. The 7? gto p is increased and 
the iterations are repeated. 

12. Obtain error term. 


21 







13. Are more answers, such as different /3, Pr, and f"(0), desired? If not, stop. 

14. The next solution is predicted by using a previously computed solution, extrapo- 
lation, or more input. If prediction of the next solution is good, the equations may be 
integrated to the largest n stop directly without going through the smaller values of 

^stop’ 
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MAIN PROGRAM FOR EXAMPLE 1 


1 

2 
3 


4 


5 

6 
7 


9 

10 

11 


12 

13 


101 

102 

103 

201 

202 

203 


DIMENSION Y( 6) »DY (6 ) » ETAEND ( 4 ) , T EST ( 4 > 

DOUBLE PRECISION Y , DY . H » ET A » X . B 
EXTERNAL DIFF 
COMMON B 

READ (5.101) H.X.B.DELPR.ETEST 

READ (5.102) N, ( ETAEND( I ) .TEST! I ) .1=1 .N) 

WRITE (6,201) H.B 
I = 1 

PRINT = DELPR 
ETA = O.ODO 

Y ( 1 ) = O.ODO 
Y(2> = O.ODO 

Y ( 3 ) = X 

Y ( 4 ) = O.ODO 
Y ( 5 ) = O.ODO 
Y ( 6 ) = l.ODO 

CALL ADAMS(6,H,ETA,0.Y,DY,1 .DIFF) 

WRITE (6,202) ETA.Y ( 1 ) ,Y ( 2 ) »Y ( 3 ) 

CALL ADAMS (6,H,ETA,1,Y»DY»1»DIFF) 

IFtETA.LT. PRINT) GO TO 6 

PRINT = PRINT + DELPR 

WRITE (6,202) ET A , Y ( 1 ) , Y ( 2 ) » Y ( 3 ) 

I F ( ETA.LT.ETAENDf I ) ) GO TO 4 

DELX = (Y(5)*(1.0D0— Y(2) ) — Y ( 6 > *Y ( 3 > ) / ( Y < 5 ) **2+Y ( 6 ) **2 ) 
E = ( 1 .ODO-Y (2 ) 1**2 + Y ( 3 ) **2 
X = X + DELX 

IF( ABS(DELX/X) .GT.TESTf I ) ) GO TO 3 
IF< E.LT.ETEST ) GO TO 12 
IF(I.EQ.N) STOP 
I = 1 + 1 
GO TO 3 

WRITE (6,203) E 
READ (5,103) B 
WRITE (6,201) H.B 
GO TO 3 

FORMAT ( 2D10.0,D20.0»2E10.0 ) 

FORMAT ( 12,/, (2E10.0) ) 

FORMAT (D20.0) 

FORMAT ( 3H1H= . D1 6 . 8 , 5X , 2 HB= , D20 . 10) 

FORMAT (4D25.8) 

FORMAT ( 3H E=»E16.8) 

END 



C SUBROUTINE TO COMPUTE THE DIFFERENTIAL EQUATIONS OF EXAMPLE 1 

SUBROUTINE D I FF { ET A , V , DV ) 

DIMENSION Y ( 1 ) »DY< 1 ) 

DOUBLE PRECISION ETA * Y *DY *B 
COMMON B 
DY( 1) = Y ( 2 ) 

DY ( 2 ) = Y ( 3 ) 

DY ( 3 ) = — Y ( 1 ) *Y ( 3 ) + B*( Y<2 >**2-1.0D0) 

DY ( 4 ) = Y ( 5 ) 

DY( 5) = Y ( 6 ) 

DY ( 6 ) = — (Y<4)*Y(3)+Y( 1 )*Y(6) ) + 2 . ODO*B*Y ( 2 ) *Y ( 5 ) 

RETURN 

END 
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MAIN PROGRAM FOR EXAMPLE 2 


1 

2 
3 


4 


5 

6 
7 


9 

10 

11 


12 

13 


101 

102 

103 

201 

202 

203 


DIMENSION Y ( 1 5 ) , DY ( 1 5 ) , ETAEND ( 4 ) , T EST < 4 ) 

DOUBLE PRECISION Y , DY , H » ETA , X ,Z 
EXTERNAL DIFF 
COMMON PR 

READ (5.101) H.X.Z.PR.DELPR.ETEST 

READ (5.102) N, (ETAENDt I ) .TEST t I > , 1=1 ,N> 

WRITE (6.201) H.PR 
I = 1 

PRINT = DELPR 
ETA = O.ODO 
Y ( 1 ) = O.ODO 
Y ( 2 ) = O.ODO 

Y ( 3 ) = X 

Y ( 4 ) = l.ODO 

Y ( 5 ) = Z 

Y ( 6 ) = O.ODO 
Y ( 7 ) = O.ODO 
Y ( 8 ) = l.ODO 
Y ( 9 ) = O.ODO 
Y ( 10 ) = O.ODO 
Y ( 1 1 ) = O.ODO 
Y( 12) = O.ODO 
Y ( 1 3 ) = O.ODO 
Y ( 14 ) * O.ODO 
Y ( 1 5 ) = l.ODO 

CALL ADAMS (15.H.ETA.0.Y.DY.1. DIFF) 

WRITE (6.202) ETA. Y ( 1 ) , Y ( 2 ) «Y ( 3 ) . Y ( 4 ) * Y ( 5 ) 

CALL ADAMS ( 1 5 . H , ETA , 1 , Y . DY . 1 , D I FF ) 

I F ( ETA. LT. PRINT) GO TO 6 
PRINT = PRINT + DELPR 

WRITE (6,202) ETA»Y(1),Y(2),Y(3)»Y(4),Y(5) 

I F < ETA.LT.ETAENDf I ) ) GO TO 4 

All = Y(7)**2 + Y ( 9 ) **2 + Y(8)**2 + Y(10)**2 

A12 = Y(7)*Y(12) + Y ( 9 ) *Y ( 14 ) + Y(8)*Y(13> + Y(10)*Y(15) 

A21 = A12 

A22 = Y ( 1 2 ) **2 + Y ( 1 4 ) ** 2 + Y(13)**2 + Y(15>**2 
B 1 = - ( Y ( 2 ) * Y ( 7 ) + Y ( 4 ) * Y < 9 ) + Y ( 3 ) *Y ( 8 ) + Y(5)*Y<10>) 

B2 = — ( Y ( 2 ) * Y ( 1 2 ) + Y ( 4 ) * Y ( 14 ) + Y(3)*Y(13) + Y(5)*Y(15)) 
DEM = All *A2 2 - A21*A12 
DELX = (B1*A22-B2*A12)/DEM 
DEL2 = (B2*A11-B1*A21 ) /DEM 

E = Y(2)**2 + Y ( 3 ) **2 + Y(4)**2 + Y(5)**2 
X = X + DELX 
Z = Z + DELZ 

I F ( ABS ( DELX/X ) .GT .TEST ( I ) .QR.ABSI DELZ/Z ) .GT .TEST ( I ) ) GO TO 
IF(E.LT.ETEST) GO TO 12 
IF(I.EQ.N) STOP 
1 = 1 + 1 
GO TO 3 

WRITE (6,203) E 
READ (5,103) PR 
WRITE (6,201) H.PR 
GO TO 3 

FORMAT ( 3D10.0.3E10.0 ) 

FORMAT ( 12,/, (2E10.0) ) 

FORMAT (E10.0) 

FORMAT (3H1H=,D16.8,5X»3HPR=,E16.7) 

FORMAT (6D18.8) 

FORMAT ( 3H E=,E16.8> 

END 



SUBROUTINE TO COMPUTE THE DIFFERENTIAL EQUATIONS OF EXAMPLE 2 
SUBROUTINE D I FF ( ET A , Y . DY ) 

DOUBLE PRECISION ETA»Y»DY 
DIMENSION Yt 1 ) *DY ( 1 ) 

COMMON PR 
DY I 1 ) = Y ( 2 ) 

DYI 2) = Y ( 3 ) 

DY ( 3 ) = 2 »ODO*Y (2)**2-3.OD0*Y(l)*Y(3)-Y(4) 

DYI A) = Y ( 5 ) 

DY I 5 ) = -3.0D0*PR*Y( 1)*Y(5) 

DYI 6) = Y I 7 ) 

DYI 7) = Y I 8 ) 

DYI8) = 4.0D0*Y(2)*Y(7)-3.0D0*(Y(3)*Y(6)+Y<l)*Y<8n— YI9) 

DYI 9) = Y I 10 ) 

DYI 10) = -3.0D0*PP*IY(5)*Y(6 )+Y I 1 )*Y I 10) ) 

DYI 11) = Y I 1 2 ) 

DYI 12) = Y< 13) 

DY I 13) = A.OD0*Y<2)*Y(l2)-3.0D0*( Y I 3 ) *Y I 1 1 )+Y 1 1 ) *Y 1 1 3 ) )-Y(lA) 
DYI 14) = Y I 1 5 ) 

DYI 15) = -3.0D0*PR*( Y( 5)*Y( 11 ) +Y I 1 )*YI 15) ) 

RETURN 

END 
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MAIN PROGRAM FOR EQUATIONS CONTAINING A PARAMETER 


1 

2 
3 


4 


5 

6 


9 

10 

11 


12 

13 

14 


101 

102 

201 

202 

203 

204 


DIMENSION Y(9),DY(9),ETAEND(4),TEST(4) 

DOUBLE PRECISION Y , DY »H . ETA »X ,3 
EXTERNAL DIET 
COMMON B 

READ (5*101) H,X.B,BH,8END,DELPR.ETEST 
READ (5.102) N, ( ETAENDl I ) »TEST( I ) . 1=1 *N) 

WRITE (6,201) H 
I = 1 

PRINT = DELPR 
ETA = O.ODO 
Y ( 1 ) = O.ODO 
Y ( 2 ) = O.ODO 

Y ( 3 ) = X 

Y ( 4 ) = O.ODO 
Y ( 5 ) * O.ODO 
Y ( 6 ) = l.ODO 
Y ( 7 ) = O.ODO 
Y ( 8 ) = O.ODO 

Y < 9 ) = O.ODO 

CALL ADAMS(9,H,ETA,0.Y,DY,1 ,DIFF> 

WRITE (6,202) B * ET A , Y ( 1 ) , Y ( 2 ) , Y ( 3 ) 

CALL ADAMS (9 *H,ETA,1,Y,DY,1,DIFF) 

IF( ETA.LT. PRINT) GO TO 6 

PRINT = PRINT + DELPR 

WRITE (6,203 ) ET A , Y ( 1 ) , Y ( 2 ) , Y ( 3 ) 

I F ( ETA.LT. ETAEND ( I ) ) GO TO 4 

DEM = Y(5)**2 + Y ( 6 ) **2 + Y(8)**2 + Y(9)**2 

DELB = -( (Y(2)— l.ODO) *Y ( 9 )+Y (3)*Y(9))/DEM 

DELX = -( ( Y ( 2 )-l .ODO ) *Y < 5 )+Y ( 3 )*Y ( 6 ) ) /DEM 

E = ( 1 .ODO-Y ( 2 ) ) **2 + Y ( 3 ) **2 

X = X + DELX 

B = B + DELB 

IF ( ABS( DELX/X) .GT .TEST ( I ) .OR.ABSI DELB/B) ,GT .TEST ( I ) ) 
I F ( E.LT.ETEST) GO TO 12 
I F < I.EQ.N) STOP 
1 = 1 + 1 
GO TO 3 

WRITE (6,204) E 
I F ( B.GT.BEND) STOP 
B = B + BH 

X = X - BH* ( Y ( 5 ) *Y ( 8 ) +Y ( 6 ) *Y< 9 ) ) / ( Y ( 5 ) **2 + Y (6 ) **2 ) 

GO TO 3 

FORMAT (3D10.0.4E10.0) 

FORMAT ( 12,/, (2E10.0) ) 

FORMAT (3H1H=,D16.8) 

FORMAT ( 3H B= , D30 . 1 6 , / , 4D2 5 . 8 ) 

FORMAT (4D25.8) 

FORMAT ( 3H E=,E16.8) 

END 


GO TO 3 
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C SUBROUTINE TO COMPUTE THE DIFFERENTIAL EQUATIONS OF, EQUATIONS 

C CONTAINING A PARAMETER, EXAMPLE 

SUBROUTINE D I FF ( ETA , Y , DY ) 

DOUBLE PRECISION ETA * Y*DY »B 
DIMENSION Y(l) ,DY(1) 

COMMON B 
DY ( 1 ) = Y ( 2 ) 

DY ( 2 ) = Y ( 3 ) 

DY ( 3 ) = B*(Y(2 )**2-1.0) - Y ( 1 ) *Y ( 3 ) 

DY ( 4 ) = Y ( 5 J 
DY ( 5 ) = Y t 6 ) 

DY ( 6 ) = 2.0*B*Y(2)*Y(5)-Y(4)*Y(3)-YU)*Y<6) 

DY ( 7 ) = Y ( 8 ) 

DY ( 8 ) = Y ( 9 ) 

DY ( 9 ) = -<Yd)*Y(9)+Y(7)*Y<3> ) + 2 . 0*B*Y ( 2 ) *Y ( 8 ) + Y<2)**2-1.0 

RETURN 

END 
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APPENDIX C 


INTEGRATION SUBROUTINE 
Description of Flow Chart 

The numbers in the following description refer to flow chart box numbers (fig. 8) and 
formula numbers of FORTRAN listing. 

1. Test to determine if subroutine is to be initialized or if it is to advance the solu- 
tion. 

2. Test if subroutine is to use Runge-Kutta or Adams-Moulton integration scheme. 
Subroutine always uses Runge-Kutta for the first three increments. 

3. Counter is initialized for Adams-Moulton scheme. 

4. Counter is initialized for Runge-Kutta scheme. 

5. Derivatives are evaluated at the initial value of the independent variable. 

6. Counter is tested. 

7. Counter is incremented by one. 

8. Counter is incremented by one. 

9. Counter is incremented by one. 



Figure 8. - Flow chart of integration subroutine. 
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10. Derivatives are evaluated at the proper values of the independent variable, and 
Runge-Kutta formulas are applied. 

11. Derivatives are evaluated, and Adams-Moulton formulas are applied. Correc- 
tion equation is used once for each increment. 

12. Solution is advanced one increment. 

This subroutine requires that a subroutine to compute the differential equations be 
supplied. A discussion of integration equations used in this subroutine is given sub- 
sequently. 
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non 


Integration Subroutine 


INTEGRATION SUBROUTINE FOR INTEGRATING SYSTEMS OF DIFFERENTIAL 
EQUATIONS OF THE TYPE DY = F(X,Y). BY RUNGE-KUTTA OR ADAMS-MOULTON 
FORMULAS 

SUBROUTINE ADAMS ( N , H , X , I S ET . Y , DY ♦ I NDEX * F ) 

DOUBLE PRECISION P , Y , YR » DYLL , DYLL L , DYL , DY » DYR , C2 ,C3 » CA , X , H 
DIMENSION Y( 1) »DY ( 1 ) *P ( 24 J » YR < 2A ) »DYLLL ( 2 A ) *DYLL ( 2A ) ,DYL t 2A) » DYR ( 
*2 A ) , C2 ( 2 A ) »C3(2A) » CA ( 2 A ) 

1 I F ( ISET.GT.O) GO TO 6 

2 I F ( INDEX. LE.O) GO TO A 

3 K = 2 
GO TO 5 

A K = 1 

5 CALL F ( X * Y » DY ) 

RETURN 

6 GO TO (10.7,8.9.11) ,K - 

7 K = 3 

GO TO 10 

8 K = A 

GO TO 10 

9 K = 5 

10 DO 1001 1=1, N 

1001 PCI) = Y ( I ) + ( H/2 • ODO ) *DY ( I ) 

CALL F (X + H/2.0D0.P.C2 ) 

DO 1002 1=1, N 

1002 PCI) = Y C I ) + (H/2.0D0)*C2< I ) 

CALL F (X+H/2.0DO.P.C3 ) 

DO 1003 1 = 1, N 

1003 PCI) = Y ( I ) + H*C3 ( I ) 

CALL F(X+H,P»CA) 

DO 100A 1 = 1, N 

1 0 0 A YR ( I ) = Y ( I ) + ( H/6.0D0 ) *< DY ( I ) +2 ,0D0*C2 ( I )+2 .0D0*C3 ( I ) +CA ( I ) ) 

CALL F ( X+H , YR , DYR ) 

GO TO 12 

11 DO 1101 1 = 1, N 

1101 PCI) = Y ( I ) + ( H/2A.0D0 ) * ( 55,0D0*DY ( I ) -59. ODO+DYL ( I ) +3 7 . ODO+DYL L 

* ( I )-9.0D0*DYLLL ( I ) ) 

CALL F ( X + H , P , DYR ) 

DO 1102 1 = 1, N 

1102 YR C I ) = Y ( I ) + <H/2A.0D0)*(9.0D0*DYR( I >+19.0D0*DY( I )-5.0D0*DYL( I ) 

* +DYLL ( I ) ) 

CALL F (X + H.YR .DYR ) 

12 X = X+H 

DO 1201 1=1, N 

Y ( I ) = YR ( I ) 

DYLLL ( I ) = D YLL ( I ) 

DYLL ( I ) = DYL ( I ) 

DYL ( I ) = DY < I ) 

1201 DY ( I ) = DYR ( I ) 

RETURN 

END 
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Integration Formulas 


The formulas described here are used in the integration subroutine given in appen- 
dix B. Let the system of n first-order differential equations to be solved be given in 
the form 


y* = f(x, y) 

where f, y, and y' are n dimensional vectors. With initial conditions 

y(x G ) = y 0 

let h be the increment of the independent variable x, with the notation Xj = x Q + jh. 
Then the Runge-Kutta formulas are (ref. 10): 

yj+i = yj + + 2E i + 2 h + 

where k^, kj, kg, and kg are n dimensional vectors and 

k Q = hf(xj,yj) 

E i = “( x i+V2’yj + i E o) 

E 2 = “( x i + i/2’yj + l 5 i) 

and 

k 3 = hf(x j+ 1 ,y j + k 2 ) 

The Adams-Moulton predictor -corrector formulas of reference 10 are given by 

= + S ( 55i S ' 59? i-! + 37? i- 2 ~ ^i- 3 ^ 

= + Yi + 19? i ‘ 5? J-» + ? i- 2 ) 
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where 


y j+l = ^ x j+l ,y ]>l) 

and the superscripts P and C stand for predicted and corrected values, respectively. 
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